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SUMMARY 

In this work, the three dimensional temperature distribution 
in the cooled rotor of a radial inflow turbine is determined 
numerically using the finite element method. Through this 
approach, the complicated geometries of the hot rotor and 
coolant passage surfaces are handled easily, and the temperatures 
are determined without loss of accuracy at these convective 
boundaries. Different cooling techniques with given coolant 
to primary flow ratios are investigated, and the corresponding 
rotor temperature fields are presented for comparison. The data 
obtained from the present analysis were found to be in agree- 
ment with the available experimental measurements. 

The present work can easily be used in combination with a 
finite element stress analysis, to investigate the thermal 
stresses corresponding to the different cooling arrangements. 

This can provide valuable information concerning the critical 
locations of possible creep, rupture or fatigue, for a given 
centrifugal, thermal and aerodynamic loading. 






INTRODUCTION 


The radial inflow turbines offer many advantages over axial 
flow turbines in small gas turbine applications. Besides 
operating at higher efficiency and yielding greater temperature 
drop and pressure ratio per stage, the radial inflow turbine 
rotor can be cast at a relatively low cost. Further improvements 
in gas turbine engine efficiencies require increased gas 
temperatures at the turbine inlet. As a result of metallurgical 
limitations, higher gas stream temperatures can be permitted 
without reducing the allowable stress levels through effective 
cooling of the turbine rotor. 

The various aspects of axial flow turbine cooling have 
been thoroughly investigated and the different techniques 
adequately developed. Several experimental and theoretical 
studies dealing with axial flow turbine cooling can be found 
in the literature. The merits and limitations of the different 
cooling techniques, namely convection, transpiration, and 
film cooling are discussed in Reference [1] . Most internal 
cooling systems have been designed using semi-empirical methods 
to achieve the highest possible effectiveness. Reference [2] 
presents a review of the present state of the art for the 
internal cooling of turbine nozzles in aircraft applications. 

Although radial turbine development, has progressed up to 
the limits of stress operating conditions, its rotor cooling 
did not achieve the high level of sophistication accomplished 
in the axial turbines. As a matter of fact, until recently 
very little research work has been reported dealing with 
radial inflow turbine rotors. Branger [3] investigated 
experimentally the effectiveness of veil cooling the hub side 
of the rotor. Ee found that the cooling effectiveness was 
larger at the rotor tip, and decreased as the cooling film is 
heated and mixed with the hot turbine flow. Petrick and 
Smith [4] measured the temperatures of a radial inflow turbine 
rotor which was cooled from its backside. While veil and 
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rotor backside cooling techniques are effective for the rotor 
disk, they induce little variation in the blade temperature. 

Therefore, except for very highly conducting rotor materials, ••• 

only internal cooling can produce pronounced reduction in the 
blade temperatures as reported in Reference [51. 

Due to the high cost and difficulties associated with 
experimental studies , computational methods can be developed 
to provide valuable information about the temperature and the 
stress distributions resulting from the different cooling 
techniques. In the present study a numerical method is developed 
to determine the temperature field .in the rotor of the radial 
inflow turbine . The computations are based on the use of the 
finite element method, with a variational statement of the 
three dimensional conduction problem. The boundary conditions 
involved at the different external and internal cooling 
surfaces of the rotor which is shown schematically in Figure 1, 
are discussed. Different cooling methods and coolant flow 
rates are investigated and the resulting cooling effectiveness 
and temperature distributions are reported. 

ANALYSIS 

The blade temperature distribution in axial flow turbines is 
often predicted from two dimensional computations at radial stations. 
Any radial variation in temperature accounts only for the spanwise 
variation in the overall convective heat transfer coefficient 
but not the spanwise conduction. In the case of radial inflow ' '! 

turbines, a similar approach cannot be used due to its complex ; 

i 

geometry. The rotor temperature distribution must therefore be j 

evaluated using three dimensional heat transfer analysis. 

In choosing a thermal analyzer for use in the case of the 
radial inflow turbine rotor , several factors such as nodal point 
placement, input efficiency , accuracy, storage requirements , 
and computer time must be taken into consideration. While the 
three dimensional finite difference thermal analyses are 
basically first order accurate, the finite element method offers 
the advantage of the capability of altering the basic accuracy 
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of the method. The finite element method was used in this 
study because of the advantages it offers in terms of input 
efficiency and nodal point placement. This is particularly 
important in the temperature computations of turbine rotors 
with internal cooling passages. Furthermore, the stress 
analysis of the turbine rotor using any of the commercially 
available finite element programs, can be greatly facilitated by 
using a thermal analyzer such as the one presented here. 

In the following, the governing equations are derived from 
variational principles. 

Governing Equations 

The field equation for the steady state three dimensional 
heat 'conduction problem in an isotropic medium can be expressed 
as : 

V* (kVT) + g = 0 Cl) 

where T is the temperature, k the coefficient of thermal 
conductivity, and g the heat generation rate per unit volume. 

The boundary conditions associated with the problem under 
consideration are : 


kVT ■ 

ii + h(T - 
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In the above equations, h denotes the convective he-at transfer 
coefficient on the surface S^, which convects heat to the flow 
at temperature T^ , and q is the specified heat flux density on 
the surface S * The union of S, and S_ forms the complete 

q b. q _ 

surface boundary S, whose outward normal unit vector is n. 

The partial differential equation (1) and its boundary 
conditions (2) and (3) can be cast in the following variational 
form according to References [6] and [7] . 


I(T) = | / [k(j§> + M§|) + k(|§) - 2g T]dV 

+ / (i h T 2 - h T w T)dS + f q T dS 04) 

S u “ 

h g 

where the desired temperature field T(x,y,z) minimizes the 
functional "I" over the domain of interest: 

5(1) = 0 (5) 

Boundary. Conditions : 

The formulation of the governing equations has been kept 
very general up till now. The surface integral was intentionally 
divided into two parts, one for the boundary surfaces over 
which a known amount of heat flux is specified and the second 
over the rest of the boundary surface convecting heat to a gas 
of known temperature. In the following section, the conditions 
associated with the different rotor boundaries will be discussed 
and the empirical expressions for the main and coolant con- 
vective heat transfer coefficients will also be given. 

Due to its rotation, any circumferent.i _, 1 nonuniformity 
in the flow conditions at inlet to the rotoi will be relatively 
averaged out, and therefore it was assumed that the inlet flow 
is axisymmetric. The rotor shown in Figure 1 was divided into 
a number of wedge sections equal to the number of blades. It 
will be sufficient to determine the temperature in any one of 
these sections since, with the assumption of axisymmetric 
inlet and exit flow conditions, the temperature field will be 
periodic. One rotor section is shown in Figure 2, with one 
rotor blade at the middle. The heat transfer was assumed to be 
negligible at the rim. With the small temperature differences 
between the blades pressure av.,d suction side, it is also! 
reasonable to assume that the heat exchange by conduction 
between two adjacent sections is negligible. Therefore, in the 
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problem formulation (Eq. 4) the heat flux q was taken as zero 
on the surface consisting namely of the rim and the two sides 
of the rotor section of Figure 2. The rest of the surfaces of 
the rotor section will be subjected to the hot gas flow. 


Main Stream Convection 

The local convective heat transfer coefficient between the 
blade and the main stream depends on the blade shape, the gas 
velocity and the position of the boundary layer transition point. 
The highest heat transfer coefficient is found at the blade 
leading edge, where the laminar boundary layer is thin, its 
value is determined from the following formula for transverse 
flow over cylinders [8] . 

Nu = 1.61 Pr^’^ Re^'^ (6) 

where Nu is the Nusselt number, Pr is the Prandtl number and 
Re is the Reynolds number based on the relative free stream 
velocity and the blade leading edge radius. 

Ainley [9] found that the mean heat transfer coefficients 
of rotating turbine blades were much higher than those for the 
same blades in two dimensional cascades. This can be attributed 
to the strong turbulence in the main flow of turbines. There- 
fore, turbulent boundary layer formulation was used in the 
present study for calculating the convective heat transfer 
coefficient between the rotor and the hot gas. This will lead 
to conservative estimates of heat transfer coefficients, resulting 

i 

in higher predicted metal temperatures, but will be used until 
more knowledge of the flow behavior in the turbine rotor becomes 
available. The following expression for the Nusselt number 
was used throughout the course of this investigations 

Pr Re (C./2) 0 * 5 

Nu = g-B (7) 

C2/C f ) ‘ +5{ (Pr-1)+S,n [1+0.8 3 (Pr-1) ] > 

Where the friction coefficient, C f , is evaluated using the 
following empirical correlation for flat plates 
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7 . C f = 0.0576 Re -0 ' 2 . 7 .( 8 ) .. 

In the case of a cooled rotor., the rest of the convective 
~ 'boundary conditions will depend on the particular cooling 

arrangement under consideration. These will be discussed in 
the following section dealing with rotor cooling. 

Rotor Cooling ? 

Because of the high centrifugal accelerations experienced 
by the rotor blades, they are cooled to lower metal temperatures 
when compared to the nozzle blades. Due to their motion 
relative to the hot gas however,- the adiabatic wall temperatV'rns 
recovered on the rotor surfaces will be lower, and furthermore 
any higher temperatures in the gases leaving the burner abov ■* 

1 the mean Values will be averaged out circumferentially 

at inlet to the rotor. Creep life consideration requires 
an accurate determination of the rotor temperature distri- 
bution. 

While Several experimental and theoretical investigations 
of turbine blade cooling can be found in the literature, very 
few investigators were concerned with radial inflow turbine 
rotor. Experimental investigations of external rotor cooling 
can be found in References [3] and [4] . in the experimental 
study of B ranger [31, the hub side of the radial rotor was 
veil cooled. It was found that veil cooling is more effective 
near the rotor tip before the cooling film is heated and mixed 
with the hot stream. Petrick and Smith- [4] measured the rotor 
temperatures when its backside is externally cooled by 
radially outward , radially inward flows and by normal impingement . 
They found that the normal impingement resulted in the highest 
cooling effectiveness while the radial inflow of coolant on the 
rotor backside gave the poorest results. 

While both rotor backside and Veil cooling are effective 
in cooling the rotor disk, they induce little variation in the 


blade temperatures for the metals normally used in rotor manu- 
facturing. Therefore, unless highly conducting rotor materials 
are used internal cooling passages in the rotor blades are 
used to obtain the desired blade temperature reduction. In 
this study the rotor temperature distribution is determined for 
the different cooling techniques shown in Figure 3. The first 
three configurations show different internal, cooling arrange- 
ments. The radial narrow holes provide the cooling passage 
in 3A. In the second arrangement 3B * which will be referred 
to thereafter as the single path, the coolant is introduced 
near the hub from the rotor backside, proceeds inside the blade 
and is discharged at the rotor tip. The cooling configuration 
C, will be referred as the double path. In this case the 

coolant is introduced at the rotor end opposite its backside, 

■** . ■ , 

cools the blade internally, turning around at the leading edge, 
^then is discharged at the blade suction side. One external 
cooling scheme is also investigated, which is shown schematically 
in Figure 3D. In the following sections, the Coolant 
temperature computations and empirical expressions used for 
the coolant convective heat transfer coefficient will be 
discussed. 


Internal Cooling : 

The flow in the internal cooling passages is affected by 
the centrifugal and the Coriolis forces. Miyazaki [10] found 
that the secondary flow is especially suppressed when the cross- 
sectional aspect ratios are further from unity. Under these 
conditions, the Nusselt number was very close to that in 
Stationary straight pipes having the same hydraulic diameter. 

The average Nusselt number for fully developed turbulent flow 
at high temperatures and heat flux densities is expressed, 
according to Reference [11] , as: 
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Nu = = 0.034 (D/L) 0,1 Pr 0,4 Re 0 * 8 (T /T ) °* 8 (9) 

is. Go 

where D is the hydraulic diameter, L the passage length, T and 
T are the coolant and hot passage temperatures respectively. 
Both Nusselt and Reynolds numbers are based on the hydraulic 
diameter, and the flow properties are referenced to T g . 

External Cooling ;: 

Patrick and Smith [4] experimentally investigated the radial 
turbine rotor backside cooling by air flowing parallel to the 
disk in both radially inward and radially outward directions 
as well as by air impinging perpendicular to the disk. They 
found that the radially outward flow resulted in higher values 
of convection heat transfer coefficients than the radially 
inward flow. The coolant flow perpendicular to the disk 
however, resulted in much higher values of convective heat 
transfer coefficients compared to those values obtained with 
flow arrangements parallel to the rotor backside. They also 
used their experimental measurements to derive empirical 
expressions for the convective heat transfer coefficient in 
all -three cases. The experimental data of this reference is so 
scattered that the authors themselves recommended limiting 
it application only to the range of the parameters in their 
study. 

Experimentally determined average Nusselt number for 
radially outward flow on a rotating shrouded disk are reported 
by Kaynes and Owen [12] . The influence of the coolant flow 
rate and the shroud and backside clearances were found to be 
less pronounced at high rotational Reynolds numbers. The 
Nusselt number would approach the free rotating disk Values 
for the high speeds and low coolant flow rates which are 
involved in radial turbine applications. , The free disk 
correlation was therefore used in the present study to deter- 
mine the local Nusselt number variation along the rotor backside 
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with the cooling air flowing radially outward. Kreith et al. 
[13] found that the critical Reynolds number of enclosed 
rotating disks with source flows is lower than the free disk 
values. Therefore, neglecting any small laminar core, that 
can exist on the rotor backside, the local Nusselt number at 
any radius r, was calculated using the following equation [14] 


Nu 



2 0.8 

= 0.0195 (^-) 


( 10 ) 


where ta is the angular velocity of the rotor, and v is the 
kinematic viscosity. 


Coolant Temperature Computations : 

The temperature of the cooling stream increases along its 
flow path due to the heat exchange with the hot rotor. This 
in turn affects the cooled surface temperature as well as the 
convective heat transfer coefficient. The computations of the 
rotor and coolant flow temperature distribution must therefore 
be performed concurrently. The rotor temperature computations 
are carried out using the finite element formulation of the 
variational statement for the three dimensional conduction 
problem. The cooling air temperature is evaluated in a 
separate program in order to achieve the desired accuracy 
without significantly increasing the storage requirements of 
the conduction problem. 

A simple energy balance equation was used to determine 
the coolant temperature rise, AT , over an incremental length, 
AL, of the passage. 


AT „ == | — h(T - T ) AL 

c mC s c 
P 


( 11 ) 


where P is the perimeter of the cooling passage, and m, the 
coolant flow rate. 


The same procedure was followed in the rotor backside 
coolant stream, to determine the variation in its temperature 
over an increment of radius, Ar. 

AT = h(T - T ) Ar (12) 

c me s c 
P 

The convective heat transfer coefficients in the above equations 
are evaluated by using equations (9) and (10) respectively. 

Finite Element Representation : 

The solution domain is discretized into a number of three 

dimensional finite elements. The temperature variation within 

(e) 

each element T (x,y,z) , is generally represented by an 
equation of the form [15] : 

T (e) (x,y,z) = {N} fc {T} <e} (13) 

4 » 

Where {N} is the row vector of the interpolation functions 

/ 0 ^ 

which depend on the nodal coordinates, and {T} is the 
column vector consisting of the nodal values of the temperature 
associated wiih that element. 

The interpolation functions N^, are chosen to satisfy 
continuity requirements to ensure the convergence of the 
solution. The integral over the whole domain in equation (4) 
can therefore be represented as the sum of the integrals over 
all the elements: 

e=M , . 

I = I I (e) (14) 

e=l 

Where M is the total number of elements in the solution domain. 

It is required to stationalize I (T) with respect to 
all the nodal values of T in the solution domain, however from 
equation (14) we can write 



. -i .. -.*Mr — H ... 
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51 


0 


(15) 


(o) ( 

I is taken only with respect to the 

with the individual element (e) . 


1, 2, ... r (16) 

number of nodal temperatures per element. 
It would therefore be necessary to derive only the single 
element equation. 

The simple four node tetrahedral elements , with linear 
interpolation functions were used in this analysis. This 
choice was dictated by computer storage considerations. In 
spite of their simplicity, these elements are very versatile 
however because of the ease with which they can be assembled 
to fit the complex three dimensional geometries with reasonable 
fidelity. A comparison of these elements with other higher 
representations 116], shows the improvement in accuracy 
is not great. The derivation of the single element equation 
for the particular finite element and interpolation function 
used in this study is given in Appendix A. When these 
equations for all the elements that constitute the rotor body 
are assembled, a relation of the following form is obtained: 

[K] {T} = {R} (17) 

Where [K] is the global stiffness matrix, {T> is the column 

vector of all the rotor nodal temperatures, and {R} is the 

column vector of thermal loading. The finite element 

representation of the rotor section, shown in Figure 2, using 

the four node tetrahedral elements results in a set of linear 

» 

equations whose matrix of coefficients has a relatively narrow 


M 


e=l 


61 


(e) 


where the variation of 
nodal values associated 
This implies that: 


81 


(e) 


8T i 


= 0 


l = 


where r is equal to the 
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band width. It is therefore adequate to solve these equations 
directly using an algorithm based on Choleski* s decomposition 
[17], forward reduction and backward substitution. 

RESULTS AND DISCUSSION 

Two previous studies involving external rotor backside 
cooling and internal blade cooling that could be used to verify 
the results of this study were available in References [4] and [5] . 
Although rotor temperature measurements were reported in Reference 
[3] , the turbine geometry and flow data supplied in that 
study was not sufficient to carry out any flow or heat transfer 
computations. It was therefore decided to carry out our 
computations for a rotor and flow conditions similar to that 
given in Reference [5] except for one difference. The rotor 
disk was extended up to the blade tip to check the ability of 
our program to handle such complicated three dimensional 
geometry. 

The rotor is made of INI00, a nickel base high temperature 
alloy, and its tip diameter is 8.2 inches. The hot gas inlet 
total temperature is 2225°F at 67,000 rpm and a turbine flow 
rate of 4.9 Ib/sec. In all the cases investigated, the cooling 
air total temperature at inlet was assumed to be 850 °F. The 
resulting temperature distributions are presented in the form 
of plots of isothermal lines on the surfaces of the blade and 
its associated hub section in Figures 4 through 9. 

The temperature fields are presented for cooled as well 
as for uncooled rotors. A schematic of the four different 
cooling arrangements investigated is shown in Figure 3. The 
computations were executed on an IBM 370 time sharing system. 

The computation time depended on the number of nodal points 
and on the particular internal cooling passage investigated , 
which in turn affects the band width of the stiffness matrix# 

For the simpler cases of solid rotor blades, the computer time 
was 15 seconds. Aside from the uncooled rotor, this includes 
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the cases of external rotor backside cooling, and the first 
internal cooling arrangement. With the very narrow cooling 
passages of that cooling arrangement, heat sinks were 
introduced to represent the amount of heat absorbed by the 
coolant in the elements which included any portions of the 
passages. The discretization of these cases involved 256 
nodal points, 675 finite elements. The number of nodal 
points and elements were larger for the rotor with internal 
cooling passage, in which the discretization involved the 
modeling of the coolant passage. The case involving rotor 
cooling through the single path shown in Figure 3 was 
modeled using 300 nodal points, and 837 elements. Larger 
number of nodes and elements, 380 and 1038 respectively, 
were used in the more complex double path cooling resulting 
in larger stiffness matrix band width and longer computational 
time. Additional external convective surface elements and 
larger stiffness matrix band widths, were naturally involved 
in the last two cases. 

The temperature field in an uncooled rotor is shown in 
Figure 4. The figure shows two views of the rotor section 
as seen from the blade pressure and suction sides. As 
expected, the highest temperatures are found near the rotor 
tip, and decrease gradually towards the hub. Both centrifugal 
and aerodynamic loading cause the highest stresses at the blade 
sections near the hub. The relatively moderate radial 
temperature gradients of Figure 4 are not expected to contribute 
significantly to the stress field. Rotor cooling should be 
expected to reduce the high temperatures of 1550 °F near the 
blade hub. If the radial temperature gradient resulting from 
a particular cooling arrangement is considerably high, it can 
augment the stresses produced by centrifugal and aerodynamic 
loading. These two factors have to be taken into consider ation 
when the different cooling configurations are compared. The 
losses incurred by the coolant injection after circulating 
through the particular cooling path is another important factor 
to be considered [18] . This however is beyond the scope of 
this study. 
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Figure 5 shows the isothermal lines for 1.5% internal 
cooling through the single path in the rotor blades. It is 
clear that even for this small amount of coolant mass flow, 
a reduction in the temperatures of 100 to 150 °F is achieved 
in the highly stressed blade regions near the hub. The 
blade temperatures were computed for the same cooling 
arrangement in Reference 15] , and the resulting temperature 
fields were given in Figures 81 and 82. When the blade 
temperatures of Figure 5 are compared with the computed results 
of Reference [5] , lower blade temperatures can be observed 
in the latter data. We must point out however that, the 
blade temperatures of Reference [5] were first iterations, 
calculated assuming zero heat flux at the blade hub. It was 
mentioned in that reference that the adiabatic blade end wall 
can account for about 50°F in lower calculated temperatures 
at the smaller radius hub sections. 

The temperature field for 3% internal blade cooling 
through a double path is shown in Figure 6. In this arrange- 
ment, 0.5% cooling is discharged at the tip to avoid choking. 
The computed rotor temperature field for this cooling arrange- 
ment is shown in Figure 6, When the blade temperatures of this 
figure were compared with those in Figures 89 and 90 of 
Reference [5] , it was found that they are in close agreement. 
The results of the last reference for this cooling arrangement 
were obtained after three iterations, to account for the 
rotor blade end wall heating effect at the hub. It can be 
seen from Figure 6 that the double path cooling arrangement 
results in a considerable temperature reduction in the blades 
highly stressed areas. This temperature reduction is achieved 
however with large temperature gradients in the blade that 
can produce considerable thermal stresses. 

In the two internal cooling arrangements just described, 
the lowest metal temperatures can be observed around the area 
of the coolant introduction in the rotor section. That can 

J 

explain the considerable reduction in the blade highly stressed 
regions in the case of a single path cooling even at the 
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relatively low coolant mass flow of 1.5%. With 3% double 
path cooling, the temperature reduction was not great, since 
the coolant is introduced at the opposite side of the rotor, 
going through a longer path before reaching these regions. 

In the third internal cooling arrangement investigated, 
the coolant path consists of a number of holes drilled radially 
in the rotor blade. The resulting temperature distribution 
with 1.5% coolant is shown in Figure 7. It can be seen that 
this arrangement with five cylindrical cooling passages of 
0.06 inch diameter results in the maximum overall reduction 
in the metal temperature, as well as in the lowest blade 
temperature near the highly stressed regions. Although this 
is desirable for better creep life, it is obvious that large 
radial temperature gradients prevail along the rotor section. 
Furthermore, large stress concentrations around the narrow 
radial cooling passages are associated with this arrangement. 

Additional results are shown in Figures 8 and 9, which 
were obtained for 1.5% and 3% external rotor disk backside 
cooling by radial outflow. It can be observed that the rotor 
backside temperatures are almost invarient near the axis and 
up to a radial distance greater than half the tip diameter, 
then increases sharply towards the tip . The temperature 
distribution of Figure 9 can be qualitatively compared with 
the experimental data given in Reference [4] for the corresponding 
cooling arrangement. Although the ratios of coolant to main 
flow investigated experimentally in Reference [4] were generally 
high, the temperature measurements at the lowest coolant mass 
flow of 4% showed the same tendencies as our temperature field 
computations . 

The cooling effectiveness was computed with 1.5% and 3% 
cooling mass flow for the rotor disk external cooling arrange- 
ment. The effectiveness, n, was defined according to the 
following expression: 
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Where T rc and T rh are the rotor temperatures with and Without 
cooling, respectively. It can be seen from Figures 10 and 11 
that, as expected, the maximum effectiveness occurs at the 
point of impingement and decreases towards the rotor tip. 

By comparing Figures 10 and 11 with Figures 8 and 9, it is 
seen that the effectiveness- curves are very similar to the 
isothermal curves. Therefore, it was unnecessary to present 
additional figures showing the effectiveness for the other 
cooling arrangements , since the inlet coolant temperature , 

T , was kept the same in all the cases considered. 

If Figures 8, 9 and 4 are compared, it can be seen that 
the rotor backside cooling is more effective in cooling the 
rotor disk, but leaves the blade temperatures practically 
unaffected. Therefore, this cooling arrangement cannot be 
effective in reducing the blade temperatures except for highly 
conductive rotor materials. It is clear however that it is 
advantageous from the stresses point of view since it mainly 
results in axial temperature gradients. Using rotor backside 
cooling in combination with internal cooling might therefore 
offer some advantages . This can be especially true if the 
internal cooling air is introduced at the rotor backside, as in 
the case of the single path. In this case, the internal cooling 
air will pass initially through a rotor section With reduced 
temperatures caused by the rotor backside cooling. If half of 
a 3% coolant is used internally in a single path and the other 
half externally on the rotor backside, for example , considerably 
lower coolant temperatures can be expected in the internal 
path near the highly stressed blade regions , just as the coolant 
passage is enlarged. This combined cooling arrangement was 
investigated and the resulting temperature field is shown in 
Figure 12. It can be seen that such combination of internal 
and external cooling offer definite advantages , besides the const 
derable temperature reduction in highly stressed regions near 
the blade hub. The high radial temperature gradients in the 
blades that are present in the case of double path cooling 
with 3% coolant are avoided here. The computed temperature 
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field of Figure 12, shows mostly temperature gradients in. the 
axial direction, with the exception of the rotor disk hear 
the tip. 

Although we did not intend to present a thermal stress 
analysis in this work , some discussion of the consequences of 
the temperature fields on the rotor stress distribution was 
presented. We would like to emphasize here again that the rotor 
dimensions and flow data were taken similar to Reference [5] 
with the exception that here , the rotor disk was: extended 
radially outward up to the tip. Although this will not 
affect the temperature field, it would naturally, result in a 
different stress distribution than in Reference [51 . Our 
computations showed that while internal cooling results in 
lower blade temperatures , particularly at the highly stressed 
regions near the hub, it also results in relatively high 
radial temperature gradients* This can augment the stresses 
produced by the centrifugal forces. The rotor external 
backside cooling on the other hand causes mostly axial 
temperature gradients with the exception of the rotor tip 
region where the centrifugal loading itself is insignificant. 

The combination of the two cooling schemes, namely that on 
the external rotor backside and in the internal single path, 
was found to result in the desired temperature reduction 
without the undesirable radial temperature gradient. 

The stress field produced by centrifugal, thermal and 
aerodynamic loadings can be determined using one of the finite 
element stress analysis programs. The critical locations of 
possible creep, rupture or fatigue can be determined. Based on 
such information, it can be seen whether the external rotor 
backside cooling with its axial temperature gradient is 
preferable to the internal cooling, even if double the coolant 
mass flow is needed to achieve the same metal temperature 
reduction around, the highly stressed blade regions. 


CONCLUSIONS 

A useful- numerical technique has been developed to predict 
the three dimensional temperature field in the cooled rotor of 
a radial inflow turbine. It was found that the finite element 
method is especially suitable for handling the complicated 
surface boundaries encountered in the different cooling 
arrangements. The calculated temperatures obtained using the 
present method are in good agreement with other analytical 
methods involving more tedious and time consuming computations. 
The three dimensional temperature fields, calculated using the 
present analysis were also found to agree with the available 
experimental measurements. 
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LIST OF SYMBOLS 


Skin friction coefficient- 
specific heat (Btu/lb °F) 

Internal Cooling passage hydraulic diameter (ft) 

A column vector representing the contribution of the 
convecting element to its thermal loading. Eg. (A-26) 

3 

Heat generation rate per unit volume (Btu/ft hr) 

A column vector representing the contribution of the 
heat generation within an element to its thermal 
loading, Eq. (A-16) 

Convective heat transfer coefficient (Btu/ft hr°F) 

A tensor representing the convecting elements 
stiffness matrix, Eq. (A-26) 

Coefficient of thermal conductivity (Btu/ft hr°F) 

Overall thermal stiffness matrix 

A tensor representing the conducting elements 
stiffness matrix, Eq. (A-16) 

Cooling passage length 

Total number of the finite elements in the solution 
domain 

Coolant mass flow rate (Ib/sec) 

The column vector of the elements interpolation 
functions 

Nusselt number 

Outward normal unit vector from the rotor surface 
Coolant passage perimeter (ft) 

Prandtl number 

The column vector of a surface element thermal loading 
due to a specified amount of heat flux. 
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q 

m 

r 

Re 

S 

s 

T 

{T} 

V 

V 
A 
v 


Heat, flux density (Btu/fr hr) 

Overall thermal load vector 

Radial distance along the rotor disk (ft) 

Reynolds number 

Rotor surface 

Finite element surface 

Temperature (°F) 

The column vector of nodal temperatures 

Overall volume 

Volume of finite element 

Surface area of finite element 

2 

Kinematic viscosity (ft /hr) 

Rotor angular velocity (radiance/hr) 


Subscripts 


c 

h 


3-rJ, 

k,£ 


Coolant flow 

Referring to surfaces exchanging heat with, the hot 
gas or coolant flow by convection 

Identify the various components of vectors or tensors , 
or nodes of a finite element 

q Referring to surfaces on which heat flux is specified 

including adiabatic surfaces 

s Rotor surface 

« Flow conditions in hot gas or coolant flow 

Superscripts 

(e) Refers to a finite element 

t Transpose of a tensor 
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SEC. A- A 


SEC. B-B 

DIMENSIONS IN INCHES 


FIG, 2 ONE OF THE ROTOR TWELVE UNITS 








B) RADIAL PASSAGES COOLING ' D) EXTERNAL DISC COOLING 


FIG. 3 | ROTOR COOLING CONFIGURATIONS 
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FIG, 12 ROTOR TEMPERATURE DISTRIBUTION WITH 1.5% EXTERNAL DISC COOLING 
AND 1,5% INTERNAL COOLING 



APPENDIX A 


ELEMENTS EQUATION 


In this appendix, the elements equation will be derived 
for the simple four node tetrahedral element that was used in 
this analysis. The condition to stationize, the functional I, 
will be rewritten here for an element, as given by equation (16) , 

= 0 i = 1, 2, ... r (A-l) 

When evaluating the integral over an element, it is more 
convenient to evaluate the different volume and surface 
integrals separately as follows: 


I ^ - I + i + l 

V s , s 
h q 


(A— 2 ) 


where 


1 9T^ e ^ ^ 3T^ e ^ ^ 8T^ e ^ ^ (e) 

I v = i /' [k(f§ ) + k(|i ) + k(|| ) - 2gT e ] dv 

' V * 


'3 Z 


(A-3) 


I = i / h(T^ - 2T T^ e ^)ds 
s h z A 


(A-4) 


I = / q T^ ds 

q A 


(A— 5 ) 


With the integrals in the above equations evaluated over the 
element’s volume, v, and its corresponding surface areas, A. 
With the highest temperature derivative appearing under the 
integral, being of first order, the four node tetrahedral 
elements with linear interpolation functions, which are used 
here, are the simplest elements satisfying the compatability 
and the completeness requirements. Substituting equation (A-2) 
into equation (A-l) , the elements equations can be written 
generally as: 
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(A - 6 ) 
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3T. 


31 
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3T . 
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s = 


where I , I ,1 are given by equations (A-3) to (A- 5) . 

V s h S q 

In. the following, the different volume and surface integrals 
will be evaluated separately, with the temperature variation 
within the element, given by equation (13) which is rewritten 
here as 


T (e) (x,y,z) = {N}"*" {1} 


(A- 7) 


Evaluation of the Volume Integral : 

Using equation (A-3) we can write 


= /; ik 8T ' e) 3 
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-]dv 
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1 , 2 , 3 , 4 
(A- 8) 


The derivatives in the above equation are evaluated with 
respect to each and every node associated with the element (e) . 
Generally, the four numbers assigned to the nodes of the 
tetrahedral elements are arbiterary, however for simplicity, 
and without loss of generality, the nodes are assigned the 
numbers 1 to 4 in this derivation. The linear temperature 
variation within the four node tetrahedral elements, can be 
expressed as follows [8] : 


T = N ± T ± i = 1,2, 3, 4 


(A- 9) 


With 



6v 


< a i + b i x + c i7 + a i z > 


(A— 10) 
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where 



and v is the volume of the tetrahedron defined by nodes 
k, £ in a right handed Cartesian coordinate system. 



1 

X i 

*i 

z i 

1 

V = 6 

1 

*3 


2 j 


1 

X k 


z k 


1 

X £ 


z a 


(A-U) 


i-r if 


(A-I2) 


The other constants are obtained through cycle permutation of 
the subscripts. It is possible to assume also a linear inter- 
polation function for the heat source. This is not justified 
however in our problem since the heat source term was introduced 
to represent the heat absorption in a narrow cooling passage. 
Therefore, assuming g to be constant in a particular element, 
and substituting equation (A-9) into (A- 8) we can write: 

aw t an. aw t 3N. 

'I* <!> + «>3F- 

+ k iff} 11 W - 9 N i ]av (A-13) 

The derivatives in the above integrations can be evaluated 
using equation (A-10) 
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6v 


(A-14) 


where b^, c^ and d^ are given by equation (A-ll) . 

Substituting equation (A-14) into (A-13) we can write 

31 


v 

3T. 
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= T j - 


i , j = 1,2,3 and 4 


(A-15) 


or 
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{ 3T Z} = Ek3 W " {G} 


(A-16) 


Where [k] is a four by four matrix in wiiich 
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--. ■ vj (b . b . + c • c > + d . d . ) 

36v 13 13 13 


(A-17) 


and {G} is a column vector, 


G i = ? v N i 


(A-18 } 


Evaluation of the Surface Integrals ; 

Using equation (A- 4) , we can write 


s h (e) 3T (e> 

sr “ { (h T<e) - H— 

1 A x 


(A-19) 


The linear temperature variation throughout the three node 
tetrahedral element is expressed as follows: 


T (e) = {N} fc { T} = N i T ± i = 1, 2, 3 


(A-20) 


Therefore we can write 


— ^ = / (h {N} fc {T> - h 5JN. ds 

1 A 


(A-21) 


In performing the minimisation, the convective heat transfer 
coefficient h, and the flow temperature are considered as 
invarients. Equation (A-21) , thus involves integrals of terms 
such as and , over the area of the triangular element* 

The values of such integrals are tabulated in various references 
[5, 6 and 15]. A general proof of the values of such integrals 
in one, two or three dimensional elements is given in Reference 
lA-1] . According to that reference, we can write for our 
surface elements: 


ct (3 y K I ft ! v ! 

' N 1 N 2 N 3 ds I 2A 


(A-22) 
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Where A is the area of the triangular element. Using equation 
(A-20) , we can write 


ff h N.N. ds + // h T N. ds = H,. - F. 

a 3 1 a 00 1 i: i 


i f j =1/2 and 3 (A-23) 


where 


H ij 


h A 
12 

„ H 
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i ? j 


i - 3 


(A-24) 


and F. = h T ^ 
1 «3 


i = 1, 2 and 3 (A-25) 

Using Equation (A-23) into equation (A-21) we can write 


3^ = [H] CT> - {F> (A— 26) 

i 

where [H] is a three by three square matrix, whose elements are 

given by equation (A-24) , and {F> is a column vector whose 

elements are given by equation (A-25) , The remaining term, 

3I„ /a T. , was similarly evaluated taking the heat flux to be 
h 1 

constant throughout the element, and using the linear 
temperature variation of equation (A-20.) . 


31 , 


3T. 
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(Q) — q j (I) 


(A-27) 


where 



i = 1, 2 and 3 


(A-28) 
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The surface area A of any triangle can be calculated if the 
three Cartesian coordinates of its vertices are known in a 
frame of reference. It is equal to half the magnitude of 
the vector resulting from the cross product of the vectors 
forming two sides of the triangle. The different surface 
and volume integrals over an element are given by 
equations (A-16) , (A-26) and (A-27) . If these equations 
are substituted into equations (A-l) and (A-2) . the element 
equation can be expressed as: 

Dc] {T} = (A-29) 

which involves the contribution of the integral over the volume 

of the tetrahedral element and the contribution of the proper 

surface integral over its four triangular surfaces. In 

(e) 

equation (A-29) , [k] is the four by four element stiffness 

( 0 1 

matrix, {R} is the column matrix representing the 

thermal loading, and {T} the column matrix of the nodal 
temperatures. The. thermal loading includes the contribution of 
the heat source (equation A-16) , the convecting surface flow 
(equation A-26) , and the specified heat flux q over an element 
surface (equation A-27) . 

Reference 

(A-l) Eisenberg, M.A. and Malvern, It.E. , "On Finite Element 
Integration in Natural Coordinates , 11 International 
Journal of Numerical Methods in Engineering, Vol. 7, 

No. 4, October 1975, pp. 574-575. 
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APPENDIX B 


COMPUTER PROGRAM 


The numerical solution to the resulting set of linear 
simultaneous equations is obtained, using the direct elimination 
approach in the subroutine CHOLES. It is called through 
the main program in which the computations of the stiffness 
matrix coefficients and the thermal load vector are carried 
out. This is accomplished through assembling the contribution 
of all the volume elements and the surface elements to 
equation (17) , As explained in Appendix A, all the volume 
elements contribute to the global stiffness matrix. Those 
with heat sources or sinks contribute also to the thermal load 
vector. The convective surface elements contribute to both the 
overall stiffness matrix and the thermal load vector. If the 
heat flux is specified at some surface elements, they contribute 
only to the thermal load vector. 

Although the elements equations were derived in Appendix A 
for tetrahedral elements, for the purpose of data preparation, 
the three dimensional body is discretized into pentahedral 
elements. This simplifies and reduces the size of the input 
data. Through computations in the main program each pentahedral 
element is further discretized into three tetrahedral elements 
as shown in Figure B. In the following, the program flow chart 
will be given, followed by the definition of the program 
symbols, a guide to input data preparation, then the program 
listing with sample input and output data. 
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Program Flow Chart: 


Read input data 
and echo print 


Divide each pentahedral element into 
three tetrahedral elements 


1. For each conducting element, compute stiffness 
matrix [k] and add its contribution to array A. 

2. For conducting elements with heat sources, 
compute the thermal load vector {G> and add 
its contribution to array A. 


3. For each convecting surface element, compute 
stiffness matrix [H] and add its contribution 
to array A. 

4. For each convecting surface element, compute 
the thermal load vector {F} and add its con- 
tribution to array A. 


For each surface element with specified heat 
flux, compute the thermal load vector {Q} 
and add its contribution to array A. 


Modify the array A for the nodes 
with specified temperature. 


— 

Gall Subroutine 
the simultaneous 

CHOLES to solve 
equations . 


Print output data 


r \ 

Stop 

V J 














r 


Symbol 

MTOT 

NTOT 

LTOT 

ITOT 

NSPEC 

NEW 

COND 

13, NT 

I 

ii 

I 

TRANS 

X(I) 

Y(l) 

i z<u 

t 

i * 

i KA,KB,KC 

j KD,KE,KG 

I 

S GEN 

| 

}j 

j| LA,LB,LC 


PROGRAM SYMBOLS 


Description. 

Total number of nodes in the finite element model, M. 

Total number of pentahedral elements. 

Total number of triangular elements on convective 
surfaces . 

4 

Total number of triangular elements on the body 
surface where heat flux is specified. 

Total number of nodes where the temperature is 
specified. 

Half band width of the stiffness matrix. 

Thermal conductivity of the body material, k. 

Code number, 12, can be set equal to zero or unity. 
In the first case, only the input data and the final 
temperature distribution will be printed. In the 
second, NT coefficients determined in intermediate 
computations will also be printed. 

A coordinate multiplication factor in case the 
conversion of units is necessary. 

An array of the nodes x-coordinates. 

An array of the nodes y-coordinates . 

An array of the nodes z-coordinates. 

Numbers assigned to the vertices of a typical 
pentahedral element. 

Rate of heat generated per unit volume within each 
pentahedral element. 

Numbers assigned to the vertices of a typical 
triangular surface element. 
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TINF 


The environment temperature of a surface element. 

HINF The convective heat transfer coefficient at a 

surface element. 

FLUX Specified heat flux at a surface element. 

NUM The node numbers, where the temperature is specified. 

TEMP Specified temperature value. 

T(I) An array of the temperatures at all the body nodes. 

The symbols of all program input data are explained above. 

The temporary storage variables were not defined. The one 

dimensional array A Cl) is used to store the elements of the 

stiffness matrix lower band k . . Cj ^ i ) , followed by the 

elements of the thermal load vector R. (i = 1, . . . M) . 

x 
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Preparation of Input Data : 


After the desired number of nodal points are placed through- 
out the three dimensional body and on its surfaces, they are 
designated the consecutive numbers between one and M. It is 
desirable that the discretization process simulate as closely 
as possible the original three dimensional body. The corners 
and the locations where there are abrupt changes in geometry 
or thermal boundary conditions, are obvious choices for nodal 
placement. At the regions where high temperature gradients are 
anticipated, finer discretization is employed. 

Using the nodal points as vertices, the body is then 
divided into a number of pentahedral elements. The six numbers 
assigned to the vertices of each pentahedral element are 
inputed in an order such that the first three define a 
triangular base. The following three numbers define the 
vertices of the other base taken in the same order as the 
vertices of the first base (see Figure B) . The coefficient 
of thermal conductivity and also the rate of heat generation 
(or absorption) if the element includes h-at sources (or sinks) 
are also inputs for pentahedral elements. 

The body surface is also divided into triangular elements 
using the surface nodal points. The input data of the surface 
elements depend upon the boundary conditions. The numbers of 
the nodes constituting each triangle are fed in the program 
input for each surface element. The corresponding local values 
of the film coefficient, h, and the environmental temperature, 

T , are specified for elements on convective boundaries. The 
heat flux is an input for surface elements in the region of 
specified heat flow. 

If the temperature is known at any nodal points, the numbers 
identifying such nodes and their specified temperatures are 
also given in the input. The program input format is explained 
in detail in the following section. 
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Type of Input 

Input Data 

Format 

Number of 
Cards 

Control data 

MTOT, NTOT, LTOT 
ITOT, NSPEC, NBW 

(615) 

one 

Control data 

I Z, NT 

(215) 

one 

Material properties 
St conversion factors 

COND, TRANS 

(2P10.0) 

one 

Geometrical 

X(I) 

( 8 F10.0) 

MTOT/ 8 


Y(I) 

( 8 P10.0) 

MTOT/8 


Z(I) 

( 8 F10.0) 

MTOT/8 

Topology and pro- 
perties of conduct- 
ing elements 

IA. IB. IC, ID. 
IE , IG, GEN 

(615. F10..0.) 

NTOT 

Convecting elements 
topology and 
boundary conditions 

LAZ., LBZ, LCZ , 
TINE, HINF 

(315, 2P10.0) 

LTOT 

Topology and bound- 
ary conditions of 
elements with speci- 
fied heat flux 

LAAZ , LABZ, 
LAC 2, PLUX 

(315, P10.0) 

ITOT 

Specified tempera- 
ture nodes 

NUM, TEMP 

(15, F10.0) 

NSPEC 
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C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 
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A FORTRAN COMPUTER PROGRAM 
***************** *** *** *** 


FOR SOLVING A THREE DIMENSIONAL HEAT TRANSFER PROBLEM 
****** ************ *********************************** 

USING THE FINITE ELEMENT METHOD 
******************************* 


DIMENSION X (255) ,Y (256) , Z 1 256 ) » P 1 4, 44 , S C 4* 44 ,L { 4 ) , M t.4) , Q( 4,4)., VI 4 
X,4 ) ,Wf 4,4) , AA{,4, 4 } , PC ( 44 ,G(4 ) , CC (.4) * SAC 3 ,3),H(3»3) , AC 9281 ) 

CALL UNDFLW 
IPRNT=6 

READ! 5,5 00 UMTOT, NTOT»LTCT .,t TOT, NSPEC. NBW 

READ (5 , 5 10 1 ) I Z , NT 

RE AD ( 5 , 5 20 1 ) CO ND , TRANS 

MT-C MTOT^NBW) +MTOT-C <NBW*C NBW+I i >/2) 

MTF=MT + I 
MTF^MT +MTQT 
N8Wl=NBW+l 
NB WP=N BW+2 


MTO T - 
NT □ T - 
LTOT - 
ITOT = 
NS P EC= 
NB W 

tz 


NT 

COND - 
TRANS 


NO, OF 
NO. OF 
NO, DF 
NO. OF 
NO. OF 
HALF BAND WIDTH 
CODE NUMBER FOR 


NODES 

THE BODY PENTAHEDRAL ELEMENTS 
THE SURFACE ELEMENTS OF CONVECTION 
THE SPECIFIED HEAT FLUX SURFACE ELEMENTS 
THE SPECIFIED TEMPERATURE NODES 
OF THE STIFFNESS MATRIX 
PRINTING A DESIRED NUMBER QF THE 


ARRAY 


(A) ELEMENTS AFTER EACH MAJOR STEP! WHICH IS THE CASE IF 
4 IZ* IS SET EQUAL TO UNITY) 

^DESIRED NUMBER OF THE ARRAY (A) ELEMENTS TO BE PRINTED 
=BQDY MATERIAL THERMAL CONDUCTIVITY 

-TRANSFORMATION FACTOR TO BE MULTIPLIED BY THE GIVEN 
NODES COORDINATES TO TRANSFORM THEM INTO CENTIMETERS 


THE ONE-DIMENSIONAL ARRAY <A> CONTAINS THE ELEMENTS OF THE 
STIFFNESS MATRIX (S> LOWER BAND FOLLOWED BY THE ELEMENTS 
OF THE R.H.S. VECTOR CC ) IN THE SYSTEM QF EQUATIONS CA)(T)={C) 


WRITE! 6, 6002)MTOT,NT0T, LTOT, ITOT, NSPEC, COND, NBW, IZ , TRANS, MT 


READ AND PRINT THE NODAL COORDINATES 
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30 0 0 RE AD{5 ,5002) ( X( T > , 1=1 .MTOT ) 
REtt> ( 5,5002) (Y( I ) , I=1,MTQT> 
READ(5,50U2MZ( I ) , I-1.MT0T) 

DO 1001 1=1 ,MTQT 

X( l)=TRANS*X{ I ) 

YC I )— TRANS *Y C I. ) 

Z{ I)=TRANS*Zt I ) 
tool CONTINUE 

WRITE ( 6 , 60 01 ) 

DO 2 I =1 ,MTOT 

WRITE C 6, 6003) I ,X( I ) ,Y( I > ,ZI I > 
2 CONTINUE 


ASSIGN ZERO VALUES TO ALL THE ELEMENTS OF THE ARRAY? A) 


IOC 2 DO 3 r =1 ,MTF 

a ( n=o .o 

3 CONTINUE 

IF ( NT QT« EQ • 0 ) GO TO 14 


THREE DIMENSIONAL HEAT CONDUCTION CALCULATIONS 
***** ****** *********************************** 


WRITE! 6,6101) 

I W-0 

GO TO 11 22 
100 IW = IW+1 

IFU W, EQ.NTOT) GO TO 1009 


READ AND PRINT THE NUMBERS ASSIGNED TO THE PENTAHEDRAL ELEMENT 
VERTICES AND THE RATE OF HEAT GENERATED IN IT 


1122 RE AD C 5 ,5005) IA, IB, I C, 1 D« IE,IG,GEN 

WRITE! 6 * 6600 > IA,IB,IC,ID, IE, I G, GEN 
20 0 2 IM=1 


BREAK THE PENTAHEDRAL ELEMENT UP INTO THREE TETRAHEDRAL 
ELEMENTS AND CONSIDER EACH ONE OF THEM SEPARATELY 


NKA=I A 
NK H= I D 
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NK C= IE 
NKD=IG 
GO TO 50 
30 0 NKA=IA 
NKB=IB 
NK C= I E 
NKD=IG 
GO TO 50 
40 0 NK A= I A 
NK B= I B 
NK C= I C 
NKD=IG 
C 
C 

C RE-ARRANGE THE VERTICES NUMBERS OF THE TETRAHEDRAL ELEMENT 

C AN ASCENDING ORDER 

C 
C 

50 IF(NKA-NKB) 1050,1050*1051 

1050 JK1=NKA 
JK2=NKB 

GO TO 11 50 

1051 JK1=NKB 
JK2=NK A 

115'' IF(NKC-NKD) 1G60 , 1060, 1061 

1060 JK3=NKC 
JK4=NKD 

GO TO 11 60 

1061 JK3=NKD 
JK4=NKC 

1160 IF( JKI-JK3) 1070, 1070,1071 

1070 KK 1 — J K 1 
KK2=JK2 
KK 3= JK 3 
KK 4= JK 4 

GO TO 1170 

1071 KK 1 = JK 3 
KK 2= J K 4 
KK3=JK 1 
KK 4= Ji<2 

1170 I F { KK2— KK4 >1080, 1080,1031 
1 C SO LK1=KK1 
LK2=KK2 
LK3— KK3 
LK4=KK4 
GO TO 11 SO 
1081 LK1=KK1 
LK2=KK4 
LK 3=KK 3 
LK4=KK2 

1180 IF (LK2-LK3 > 1 3 50 , 1 090, 1 190 
1090 KA=LK1 


KS=LK2 



IN 
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KC=LK3 
KD=LK4 
GO TO 2100 

1190 KA=LK1 
KB =LK3 
KC=LK2 
KD-LK4 

2 1 C 0 WR ITE { 6 , 600 4 ) K A , KB , KC * KD 
C 
C 

C CALCULATE THE TETRAHEDRAL ELEMENT CONTRIBUTION MATRIX, (K) 

C 

C 

20 0 3 P(l,l)=1.0 

P< 1,2) =X(KA ) 

P{ 1,3 )=Y (KA) 

P( 1,4)=Z(KA) 

P( 2, 1 )-l .0 
P{ 2,2)=X(KB) 

P( 2.3)=Y(KB) 

P( 2,4) — Z (KB) 

P( 3,1 ) =1 *0 
P( 3,2) =X (KC) 

P( 3, 3 )=Y (KC) 

P(3, 4)=Z(KC) 

P( 4, 1 )=1 ,0 
P( 4,2) — X (KD) 

P( 4, 3)=Y (KD) 

P( 4,4)=Z(KD) 

C0F=(P{2*2)*((P{3,3)*P(4,4>)-(P(3,4>*P<4,3})>)”(P(3,2>*(£P{2,3)*P( 
X4,4>>-(P{2»4)*P(4,3))))+(P(4,2)*({P(2,3)*P(3,4))-{P(2,4)*P(3,3)>)> 
CQG={P(1,2)*£(P{3,3)*P(4,4))-{P{3,4)*P(4,3>>))~(P(3,2)*{(P{1,3)*P( 
X4,4) ) -tP(l ,4>*P(4,3) ) ) )+ (P(4 ,2 )* ( (P( 1,3) *P{3,4) )-(P( 1 , 4) 4P ( 3,3) ) > ) 
C0H=(P(i,2)*(CP(2,3)*P(4,4))-(P(2.4)*P(4,3))))-(P(2,2)*({P(l,3>*P( 
X4,4)l-(P{l,4)*P(4,3))))+{PC4,2)tt(PC li 3) *P( 2,4) )-( PC 1,4)*P{ 2,3) ) ) ) 
COP=(P ( 1 ,2>*< (P (2,3)*P( 3 ,4) >-( P(2»4)*P<3,3 ) ) ) )-(P( 2,2) *( CP( 1 ,3>*P{ 
X3, 4) )- (P(l ,4)*P(3,3 >> ) )t(P(3,2 )*( (P( 1,3) *P { 2,4) )-( P( 1 ,4)*P(2,3 ) ) ) ) 
AD=CQF-COG+CQH-CaP 
VDL=( ABS £ AD) ) /6.0 
CT=VOL*CGND 
CTT=GEN/3.0 
CALL MINV(P,4,AF,L,M) 

S( 1 , 1 ) -0 .0 
SC 1,2) =0 ,0 
S( 1,3) -0 .0 
S(l,4)-0.0 
S( 2, 1 ) =0 .0 
S( 2, 2) = 1 ,0 
S( 2,3) =C .0 
S( 2,4) =0 .0 
S( 3, 1 )=0 «0 
S( 3,2) — 0 ,0 
S( 3,3) =1 ,0 
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5( 3*4) =0 .0 
St 4* X) -0 .0 
S( 4*2) =0 • 0 
S( 4*3) =0 .0 
St 4,4) =1 .0 

CALL MPRDC S,3, Q ,4,4 *5,0,4) 
CALL MTRAt P*V .4,4*0 ) 

CALL MPRDC V,Q.W*4,4,0 ,0*4) 
DO 1000 1=1,4 
DO 2000 J= 1*4 
AAt I * J )=CT*WC I * J) 

20 0 " CONTINUE 
10 90 CONTINUE 


C 

C 

C 

C 

r 

C 


COMPUTE THE ORDERS OF THE ARRAY (A) ELEMENTS TO 
QY ADDING THE MATRIX (K) TO THE ARRAY tA) 


QE INFLUENCED 


NQ=1 
KN=KA 
NSUM=0 
GO TO 8872 
8871 NSUM=0 

IF ( NQ • EG *2 ) KN=KB 
IFtNQ.EQ.3) KN=KC 
IF tNQ • EQ *4 ) KN=KD 
IFtNQ.EQ.S) GO TO 1929 
DO 9971 J=1,KN 
NSUM=NSUM+J 
CONTINUE 
ID G=NSUM 

IF(KN-NBWP) 7765,7766*7766 
I = IDG 

GO TO 9973 
NACC=G 

DO 7768 K=NBWP,KN 
NACC=NACC+K-N5WP+i 
CONTINUE 
1= IDG-NACC 
At I)=A ( I )+AAC NQ.NG) 

NQ=NQ+ 1 
GO TO 88 71 
ICODE- 1 
NU=KB 


8872 
997 1 

7765 
776 6 

7768 

9973 

1929 

322 5 
3226 


NV=KA 
GO TO 
NU=KC 
NV=KA 
GO TO 
NH-KC 
NV=KB 
GO TO 


3260 


3260 


3260 
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3227 NU=KD 
NV=KA 

GO TO 32 60 

3228 NU=KD 
NV=K0 

GO TO 32 60 

3229 NU=KD 
NV=KC 

3260 LNlteNU-1 
NADD=0 

DO 3261 J=1 iLNU 
NA DD=N ADD+J 

3261 CONTINUE 
ND G=NADD 

IF(NU-NBWP) 1980.1984,1981 

1980 I=NDG+NV 
GO TO 1982 

1981 JACC=0 

DO 1983 K=NBWP,LNU 
JACC=J ACC+K-NBWP-H 

1983 CONTINUE 

1 1 =NDG-J ACC 
I=I1+NV-NU+N3WP-1 
GO TO 1982 

1984 I-NDG+NV-1 

1982 GO TO (3262,3263, 3264, 3265.3266, 3267) , ICODE 


ADD THE TETRAHEDRAL ELEMENT CONTRIBUTION MATRIX (K) TO THE 
ARRAY (A) IN THE PROPER PLACES CALCULATED BEFORE 


3262 A{ I)=A(I)+AA(2,1) 

GO TO 3367 

3263 A{ I)=A ( I )+AA{ 3,1) 

GO TO 3367 

3264 A( I)=A ( I ) + AAC 3, 2) 

GO TO 3367 

3265 AC I ) ~A ( I ) 4-AAC 4,1 ) 

GO TO 3367 

3266 A ( I ) =A ( I ) + AA (4,2) 

GO TO 3367 

3267 A( I ) =A ( I )+AA( 4,3) 

3367 IC0DE= I CODE 4-1 

GO TO (9000,3225,3226,3227,3228,3229,8844) .ICODE 


C COMPUTE THE TETRAHEDRAL ELEMENT COLUMN VECTOR (G) RESULTING 

C FROM THE HEAT GENERATED WITHIN THIS ELEMENT 

C 
C 

8 844 XC = ( X( KA ) + X ( KB > +X{ KC ) +X t KD ) ) /4. 0 
YC={ Y ( KA ) -(-Y ( KB > +Y ( KC H-Y C KD ) > /4 .0 


55 


EEPRO&tJCEnUTY OF THE 
OIUEflNAIi PAGE IS POOR 


.. ...L... 


I 


ZC=(Zt KA >+Z{KB)+Z(KC)+Z{KD)>/4.0 
PC (1 ) = 1 , 0 
PC (2)=XC 
PC (3)=YC 
PC t4)=ZC 

CALL MPRDtV, PC, G, 4*4, 0.0,1) 

DO 9 I =1 ,4 
CC ( I )=CTT*Gt I ) 

9 CONTINUE 
KAP=MT+KA 
KBP=MT+KB 
KCF=MT+KC 
KDP=MT+KD 

A { KAP ) =A {K AP ) +CC{ 1 ) 

A( KBP) =A (KSP) +CC( 2 > 

At KCP ) =A (KCP ) + CC( 3 ) 

At KDP)=A t KOP ) +CC t 4 ) 

C 

C 

C RETURN TO CONSIDER ANOTHER TETRAHEDRAL ELEMENT 

C 

c 

IM=IM+ 1 

GD TO (9000 ,300,400 ,100 ) *IM 
1 0 C 9 IF(IZ.EQ.O) GO TO 14 
C 
C 

C PRINT PORTION OF THE ARRAY (A) AFTER CONSIDERING ALL THE BODY 

C ELEMENTS 

C 
C 

WRITE16* 60 OS) 

DO 12 1=1, NT 

WR ITEt 6, 6006) I * A t I ) 

12 CONTINUE 

14 IFtLTOT.EQ.O ) GO TO 23 
C 
C 
C 

C HEAT CONVECTION BOUNDARY CONDITION 

C **************** ****************** 

c 

c 

c 

WRITEt 6, 6009) 

401Q DO 15 K=1»LTQT 
C 

c 

C READ THE NUMBERS ASSIGNED TO THE CONVECTIVE TRIANGULAR ELEMENT 

C VERTICES AND THE VALUES OF BOTH THE ENVIRONMENT TEMPERATURE 

C AND THE HEAT TRANSFER COEFFICIENT 

C 
C 
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RE AD( 5 , 5006) LAZ * LBZ *LCZ • TINF,HINF 


RE— ARRANGE THE NUMBERS IN AN ASCENDING ORDER 
C 
C 

IF(LAZ-LBZ) 4012,4012,4013 

4012 LA1=LAZ 
LB 1=L3Z 
LC1=LCZ 

GO TO 40 14 

4013 LA1=LBZ 
LB 1=LAZ 
LC1=LCZ 

4014 IF(LAl-LCl) 4015,4015,4016 
40 1 5 LA 2— LA 1 

LB2=LB1 
LC2=LC 1 
GO TO 4017 

4016 LA2=LC1 
LB2=LB 1 
LC2=LA 1 

4017 IF (LB2— LC2 ) 4018,4018,4019 

4018 LA=LA2 
LB=LB2 
LC=LC2 

GO TO 40 20 

4019 LA=LA2 
LB=LC2 
LC=LB2 

40 20 WR ITE ( 6 » 60 10 1 LA, LB » LC» T INF,HINF 
C 

c 

C COMPUTE THE CONTRIBUTION MATRIX CH) OF THE CONVECTIVE ELEMENT 

C 

c 

4021 ALA={ { Y(LB)-YCLA) )*(Z{LC)-Z( LA) >-( Y(LC>-Y(LA> }*( Z(LB>-Z(LA) ) ) 
ALB=( C X(LB )-X(LA> ) *(Z( LC)-Z{LA) )-(X{LC>-X(LA)}*{Z(LB>-Z<LA>>) 
ALC={ ( X(LBJ-X( LA) > *( YCLCJ-YtLA ) )-C X( LC >- X{ LA ) ) *< Y< LB ) -Y (L A ) ) ) 
ZLA=( ABS C ALA) ) *+2,0 
ZLB=< ABS(ALB) ) **2,0 
ZLC=< ABS (ALC) >**2*0 
AREA=Q .5*1 ( ZLA + ZLB+ZLC ) * *0 .5 > 

SA (1,1 ) = 2.0 
SA<1 ,2 >=I*0 
SAU,3 )=l-0 
SA ( 2, 1 )=1.0 
SA (2*2 ) = 2*0 
SA<2,3)=1.0 
SA (3, 1 )=1*0 
5A(3,2 ) = 1 .0 
SA (3,3 > = 2.0 
CINF= ( AREA*HINF>/12*0 
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DO 16 1=1,3 { 

DO 17 J— 1*3 

H( I * J ) — C INF*SAl I , J ) ! 

17 CONTINUE 
16 CONTINUE 

C ? 

C * 

C CALCULATE THE ORDES OF THE ARRAY (AJ EL EM ENT 5 TO WHICH THE j 

C CONVECTIVE ELEMENT CONTRIBUTE i 

C j 

c ! 

NQ 1=1 ;/ 

KN1=LA ! 

NS UM=0 
GO TO 41 11 
AllO NSUM=0 

IF (NQ1 • E Q. 2 ) KN1=L8 
IF(NQ1.EQ.3) K.N1=LC 
IF (NQ1 *EQ»4 ) GO I 4400 

4111 DO 4112 J=1 *KN1 

NSUM=NSU M+J 

4112 CONTINUE 
IDGl=NSUM 

IF(KNl-NBWP) 4113*4114*41 14 

4113 I=lDGi 

GO TO 41 17 

4114 NA 0=0 

DO 4115 J=NBWP,KN1 
NAC1=NAC 1+ J-NBWP+1 

4115 CONTINUE 

1= IDG1 -NAC1 

4117 A { I ) =A ( I )+H( NO 1 * NQ 1 ) 

NQ 1=NQ 1+1 
GO TO 41 10 
44^0 I COD 1=1 
NU 1=LB 
NV 1— LA 
GO TO 4401 

4501 NU1=LC 
NV 1=LA 

GO TO 4401 

4502 NU 1 — LC 
NV 1=LB 

4401 LNU1 =NU 1—1 
NAD1=0 

DO 440 2 J=1 • LNU 1 
NAD1=NAD1+J 

4402 CONTINUE 

NO G1 =NAO 1 j 

IFtNUl-NBWP) 4403,4404,4405 

44^3 I=NDG1+NV1 I 

GO TO 4407 
440 5 JA Cl =0 
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DO 4406 J=MBWP,LNU1 
JAC1=J AC I + J-NBWP-H 

4406 CONTINUE 
J1=NDG 1-JAC1 

I = J1 + NV1-NUH-N8WP- 1 
GO TO 4407 
4.404 I=NDGl+NVl-l 
C 
C 

C ADD THE CONTRIBUTION OF THE CONVECTIVE ELEMENT TO THE PROPER 

C ELEMENTS OF THE ARRAY (A) 

C 

C 

4407 GO TO (4601,4602,4503) , ICODl 
4601 A( I)— At I )+H( 2 , 1 ) 

GO TO 4904 

460 2 A( I )=A C I > +-H{ 3 » 1 ) 

GO TO 4904 

4603 A( I)=A( I )+H(3,2) 

4904 ICCD1- ICOD l + l 

GO TO £ 90 00 ,453 1 ,4502 ,460 5) » I COD1 
C 
C 

C COMPUTE THE CONTRIBUTION VECTOR (F) OF THE CONVECTIVE ELEMENT 

C 

C 

4605 DI NF=( HI NF*T INF*AREA)/3 .0 
LAS=LA+MT 
LBS=LB+MT 
LCS=LC+MT 
C 
C 

C ADD { F) TO THE ARRAY (A) IN THE PROPER PLACES 

C 

C 

A( LAS )=A (LAS ) *DINF 
A( LBS) =A(LBS) +DI NF 
A(LCS)=A(LC5)+DINF 
15 CONTINUE 

IF(IZ.EQ.O) GO TO 23 

WRITE! 6.60 ID 

DO 21 1=1 , NT 

WRITE ( 6,6012)1 ,A( I } 

21 CONTINUE 

23 IF( ITOT.EQ.O ) GO TO 29 
C 
C 
C 

C SPECIFIED HEAT FLUX BOUNDARY CONDITION 

C ** ****** **** ************************** 

C 

C 

c 
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WR ITE( 6* 60 15) 

DO 24 K= 1 , ITQT 
C 
C 

C READ THE NUMBERS ASSIGNED TO THE SPECIFIED HEAT FLUX TRIANGULAR 

C ELEMENT AND THE SPECIFIED VALUE OF THE HEAT FLUX AT THAT 

C ELEMENT 

C 
C 

READ(5 ,5007 >LAAZ,LABZ,LACZ,FLUX 
C 

c 

C RE-ARRANGE THE VERTICES NUMBERS I N AN ASCENDING ORDER 

C 

C 

IF (LA A Z— LABZ ) 90 12, 90 12, 9013 
90 12 LAA1=LAAZ 
LABI -LABZ 
LACi=LACZ 
GO TO 90 14 
9013 LA A1=L AB Z 
LA B1=L AA Z 
LAC1=LACZ 

90 14 IF (LAA 1-LAC 1 ) 9 0 1 5 , 90 15,9016 

9015 LA A2=L AA 1 
LA82=LAB 1 
LA C2=L AC 1 
GO TO 9017 

9016 LA A2=L AC 1 
LAB2=LAB 1 
LA C2.-LAA 1 

9017 IF(LAB2-LAC2)901B, 9013, 9019 
90 18 L A A=LA A2 

LAB=LA02 
LAC=LAC2 
GO TO 90 20 
9019 LA A=LA A2 
LA E=LAC2 
LAC=LAB2 

90 20 WR ITEt 6, 60 16)LAA,LAB,L AC ,FLUX 


COMPUTE THE ELEMENT CONTRIBUTION VECTOR (Q) 


ELA=( ( Y(LAB)-Y(LAA ) )*< Z( LAC)-Z(LAA) )— (Y(LAC) — Y (LAA) }*(Z(LAB) — Z (LAA 
X ) ) )**2 =0 

BLB=( ( XtLABi-XtLAA ) ) * C Z(LAC) -ZtLAA ) >-(X< LAC) -XI LAA) ) * C Zt LA B ) -Z ( L A A 
X ) ) )**2 -0 

8LC=< ( X( LAB)-X (LAA ) )* ( Y ( LAC ) — Y (L AA ) )— ( X ( LA C ) — X ( L AA ) ) * C YC LA B ) ~Y ( L AA 
X) ) >**2 ,0 

AR EA=0 ,5 *£ ( BLA-t-BLB + BLC >**0,5) 

FI X=( AREA4FLUX )/3 *0 
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LA 1=LA A+MT 
LB I=LAB+MT 
LC1=LA C+ MT 
A( LAI ) =A (LA i ) 4-FI X 
A{ LB l ) =A {LSI) +F I X 
A( LC1 ) =A(LC1 ) +FIX 
24 CONTINUE 

IFCIZ.EQ.O) GO TO 29 
WRITE! 6*6017 ) 

DO 26 1=1, NT 

WRITE! 6, 60 18i I , A{ I > 

28 CONTINUE 

29 IF ( NSP EC ■ EQ ■ 0 ) GO TO 43 
C 

C 

c 

C SPECIFIED TEMPERATURE BOUNDARY CONDITION 

C ****************************** ********** 

C 

c 

c 

WRITE! 6, 6019 ) 

DO 30 1= l , NSPEC 
C 
C 

C READ AND PRINT THE NODE NUMBER AND ITS SPECIFIED TEMPERATURE 

C VALUE 

C 
C 

READ! 5, 5008) NUM, TEMP 
WRITE! 6, 60 20 ) NUM, TEMP 
C 
C 

C MODIFY THE ARRAY (A) BY SUBSTITUTING THE SPECIFIED TEMPERATURE 

C VALUE IN THE RESULTING SET OF SIMULTANEOUS EQUATIONS 

C 
C 

IM=NUM-1 
IP=NUM+1 
KUM1=NUM+MT 
IP I=MT + IP 

IF(NUM.EQ.l) A { NIJM ) = 1 • 0 

A( MJM1 ) = TEMP 

LSUM=0 

DO 2450 J=1,NUM 
L5UM=LSUM+ J 

2450 CONTINUE 
LDG=LSUM 

IF {NUM sEQel) GO TO 2457 
IFCNUM-NBW I) 2451,2451,2452 

2451 A(LDG>=1.0 

LD GM=LDG— NUM-Fl 
GO TO 2454 
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2452 IS UM— 0 

DO 2453 J=NBWP,NUM 
ISUM=I SUM+J-NBWP+l 

2453 CONTINUE 

LD G=UD G-ISUM 
A ( LOG ) = 1 .0 
LDGM=LDG-NBW 
GO TQ 2455 

2454 MT2=MT+IM 

DO 2456 J=MTP,MT2 
K= J— MT P+ L.D GM 
A ( J)=A (J )-( TEMP* AIK > ) 

A( K)=0 .0 

2456 CONTINUE 
GO TO 25 56 

2455 MPLC=M TP+NUM— N9W1 
MT2=MT + IM 

DO 2555 J=MPLC,MT2 
K= J— MPLC+LDGM 
A ( Ji=A ( J )-( TEMP* A { K) ) 

A ( K)=Q .0 

2555 CONTINUE 

2556 IFINUM •EQ.MTCt) GO TO 30 
JUD=MTOT-NUM 

IF (JUD .LT.NBW) GO TO 6560 
IF(NUM-NBWI) 2457.2458.2458 

2457 NUM2=NBW1-NUM 
I3=LDG 

00 646 1 K= I * NUM2 
13=1 3+NUM+K-l 
I4=NUM 1+K 

A( I4)=A( 14 )-< TEMP* A I 13 ) ) 

A ( 13) =0. 0 
6461 CONTINUE 

1 J = I 4 

IFINUM.EQ. 1) GO TO 3020 
DO 6462 K= 1 . I M 
T 3 = 1 3+ NB W 
14=1 J + K 

AI 14 ) = A I I4)-(TEMP*AtI3) ) 

A ( 13 ) =0 • 0 
5462 CONTINUE 
3020 DO 7010 J=MTP.MTF 

WR ITE ( 6*71 10 ) J .At J ) 

7110 FORMAT (4X, I3.4X.E12. 5) 

7010 CONTINUE 
GO TQ 30 

2458 I6=LDG 

DO 6463 K=1,NBW 
1 6 = 1 6+ NB W 
I7=NUM 1+K 

A ( 17 ) =A { 17) -I TEMP* A I 16) ) 

A ( I6)=0.a 


• ’ ? 
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6463 CONTINUE 
GO TO 30 
6560 I8=LDG 

DO 6561 K=1,JUD 
1 3=1 8-FNBW 
E9=NUM 1+K 

At 19 ) =A( 19) -(.TEMP* A( 13) ) 

A( ra>=o.o 

5561 CONTINUE 
30 CONTINUE 

IF(IZ.EQ.O) GO TO 43 
WR ITEt 6* 1985) 

DO 1986 1=1. NT 
WR ITEt 6*1987)1 .At I ) 

1936 CONTINUE 
C 

c 

C SOLVE THE FINAL SET OF SIMULTANEOUS EQUATIONS 

C 

C 

43 CALL CHOLEStA, MTOT, NBW1 ,0. 1. IPRNT.&9000) 

C 

c 

C PRINT THE TEMPERATURE VALUES AT ALL THE BODY NODES 

C 

C 

WRITEC 6.6023) 

DO 45 I=MTP» MTF 
J= I— MT 

WR ITEt 6, 6024) J , At I ) 

45 CONTINUE 
5001 F0RMATC6I5) 

5101 FORMAT (215) 

5201 FORMAT (2F10 .0 ) 

5QC2 FORMAT (8 FI 0.0) 

5HC5 FORMAT (6 I5.F10 .0 ) 

5006 FQ RMAT (3I5.F10.0) 

50 0 7 FORMAT t 3 I5.F 10.0) 

5003 FORMAT t I5.F10 .0 ) 

60 0 2 FORMAT { 5 5X » • MTOT=* , r 3/S5X. • NTOT=* ,I3/55X , * LTOT=* , I 3/55 X , * ITOT=* , I 
X3/55X, * NSPEC= * , I3/55X, ' C OND= » , F 5 . 2/55X , * NBW=« . I3/55X . * IZ=» . I1/55X 
X . * TRANS= * »F10 .5///45X. * S IZE OF THE STIFFNESS M ATR I X= • , 1 6 ) 

6C01 FORMATt 1H1//53X. • NODAL POINTS CO QRD I NA TE S * /30X , • NODE N0.',13X,*X 
> r . 15X * ■ Y ' » I5X, * Z*//) 

6003 FORMAT C32X. 13 »10X,F10.3,5X.F1Q.4.5X,F10*4) 

6101 FORMATt 1HI/43X ,• THREE DIMENSIONAL HEAT CONDUCTION CALCULAT IONS */4 
X 3X * * e .46( = * = ) /43X , B e , 46 { fi ** > // 61 X . * I NP UT DATA*/61X.» 

X/) 

6800 FORMAT 130X.6C4X. I 3 ) , 4X , F 1 G *4 ) 

600* FORMAT (8 1X*4(4X. 13 > ) 

60 C 5 FORMATt IH1//54X.*°0RTI ON OF THE ARRAY ( A ) ' /28X . • AFTER ADDING THE C 
XDNTRI3UTIQN OF THE CONDUCTIVE FINITE ELEMENTS*//) 

6006 FO RMA T { 54X . 15 ,4X . El 0 *3 ) 
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600 9 FORMAT 11 H1/53X.* HEAT CONVECTION BOUNOARY CON D IT I ON 0 / 49 X « * ' , 
X34!** 6 i//oiX* s INPUT DAT A ® /6i X , * « *10 ( 1 *»)// ) 

6010 FORMAT (4-5X ,3! 4X« 13), 6X * F 6 « 1 , AX * F 5 « 1 ) 

6011 FORMAT! IH1//54X, ' PORTION OF THE ARRAY { A ) • /28X , • AFTER ADDING THE C 
XONTRIBUTION OF THE CONVECTIVE FINITE ELEMENTS*//) 

6012 FORMAT ! 54X » 15 ,4X , E10 ,3 ) 

6015 FORMAT! IH1/47X, * SPECIFIED HEAT FLUX BOUNDARY CONDI TIQN* /47X , » «*3 

X8! )/47X. ' * »38( • ★' ) // 61 X, •INPUT DATA»/61X,' *, 1 C !«**)// ) 

60 16 FORMAT C5GX ,3! 4X. I3 ) ,6X,F5.1) 

60 17 FORMAT ! 1H1//54X, 'PORTION OF THE ARRAY ! A ) * /3 3X , * AFTER ADDING THE C 
XONTRIBUTION OF THE SPECIFIED HEAT FLUX ELEMENTS'//) 

60 13 FORMAT { 54X , 1 3 ,4X , E1Q ■ 3 ) 

60 19 FORMAT! 1 Hi //30X, » INPUT DATA FOR THE SPECIFIED TEMPERATURE NODES'// 
X) 

6020 FO RMAT ! SOX , 13 , 4X* F 6*1) 

1985 FORMAT (1H1//54X.* PORT! ON OF THE ARRAY t A ) ' /41X , * AFTER INTRODUCING 
XTHE SPECIFIED NODAL TEMPERATURES * // ) 

1987 FO RMAT ! 5 6X » I5» E 1 2* 5 ) 

6023 FORMAT ! 1 H1/51X, 'FINAL TEMPERATURE D I STR I BUT I ON * /SO X , * ',30 !•*»)/ 

X 50 X * * * » 30 ( * * * ) //52 X » 'NODE NUMBER* ,5X » * TEMPERATURE ' /51 X* • * , 11!'*' 
X>, 4X* * * .11! )/) 

6024 FORMAT (55X, 15, 1 IX, F12.6 ) 

9000 STOP 

END 


O 
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SUBROUTINE CHOLESC B , N, MM , I B , NT , I PRNT * * ) 

REAL*3 Sl.Tl 
DT MENS ION B( I ) 

C 
C 

c 
c 
c 
c 

c * 
c 
c 

MUD = MM-1 » 

NS = MUD*MM/2 
NM = N*MM-NS 
IF (NT-1 ) 30 ,30 ,31 
C 
c 

C BEGIN CHOLESKY ALGORITHM FOR FACTORING THE MATRIX 

C 

C 


30 

DO 

20 J = 1 ,N 


IFCJ-MUD) 1,1,2 

2 

IN 

=J— MUD 


L= 

IN J— MM) *MUD +NS 


GO 

TO 7 

1 

IN- 

= 1 


L= 

IN +IJ-l)*J/2 

7 

IF 

IJ — N+MUD1 103, 103,105 

105 

MS 

— N 


GO 

TO 10 4 

10 3 

MS 

- J+MUD 

104 

si 

= 0.0 


J1 

= J-l 


J2 

= J + l 


IF (J1 ) 4,4,3 


3 DO 5 K = I N , J 1 
T1 = B (L ) 

SI =S1 + Tl**2 


6 L=L +1 
4 T1=S(L) 

IFCT1-S1 ,LT ,0 « ) GO TO 100 
T 1 =OSQRT{ Tl - SI ) 

BID a Tl 
IF (J-N) 19,20,20 
19 DO IB 1= J 2 ,MS 
SUM = 0,0 
IF ( I —MUD )6S»68j71 
71 IN = I —MUD 

LL = IN + I-HM)*MUD +NS 
GO TO 5 
68 IN = 1 

LL = IN + t i-1 )*I /2 


{ B ) IS THE ARRAY CONTAINING THE ELEMENTS OF THE LOWER HALF BAND 
OF THE STIFFNESS MATRIX IS) 

N =NO* OF THE UNKNOWN NODAL TEMPERATURES 
MM =HALF BAND WIDTH OF THE 5TIFFNESS MATRIX +1 
SI , Tl ARE TEMPORARY DOUBLE PRSCES ION VARIABLES 
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j 


5 IFCJ1 ) 18 ,18,8 
8 IF (IN -Jl) 53, 53 ,18 
S3 DO 17 K = I N* J 1 
LM = L + K-J 

SUM = SUM + B(LL)*B(LM> 

17 LL = LL +1 

18 B(LL) = (B(LL)-SUM) /S { L ) 

20 CONTINUE 

31 NR = IB +1 
C 

c 

C BEGIN FORWARD SUBSTITUTION 


NB =NM +1 

DO 65 K =1 *NR 

B ( NB ) =B(NB)/3{1> 

DO 60 I = 2,N 
IF ( I -MUD ) 21,21,22 

22 IN ~ I -MM 

KS = IN* MUD +NS 
M5 =MUD 
GO TO 27 
21 IN = 0 
MS = 1-1 
KS = M5*I /2 
27 SUM = 0.0 

DO 61 J— 1,M5 
JR = J+IN 
L = JR +KS 
JR = JR + NB - 1 
61 SUM = SUM +B(L) *B{JR) 

ID = I + KS 

60 B ( JR + 1 ) = (B(JRF1) —SUM ) /B (ID) 
65 NB = NB +N 
C 
C 

C BEGIN BACKWARD SUBSTITUTION 

C 

c 

NO ■- N M + N 
DO 75 K = 1 , NR 



S{ 

NB > 

— 

B( 

NB) 

/B(NM ) 


DO 

30 

II 

— 

2 , N 



I 

= N 

-I 

r 

+ 1 



IF 

{ I- 

MUD 

)4 1 , 4 1 

,95 

95 

ID 

= 1 

+ ( 

I- 

MM)* 

MUD +NS 


GO 

TO 

42 




At 

ID 

= 1 

+ ( 

I- 

1 )*I 

/2 

42 

IF 

(I- 

N+MM) 

43, 

43,45 

45 

MS 

= 1 

1-1 





GO 

TO 

76 




A3 

MS 

- 

MUD 





■Li 
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76 SUM =0.0 

DO 81 J =L ,M5 
JR = I +J 

rF (JR— MUD) 98,98,99 
99 l = I +-( JR— MM ) *MUD + NS 
GO TO 82 

98 L = I +{ JR— 1 ) * JR/2 

82 JR=NB-N+JR 

81 SUV =SUM + B(L)*B(JR) 

JR =NB -n + 1 

80 B{ JR ) = (S(JR)-SUM) /B(ID) 

75 NB =NB + N 
RFTURN 

1^0 WRITE! I PRNT ,93 1 ) 

901 FORMAT (1H1, 1 THE MATREX IS NOT POSITIVE DEFINITE IN THIS PROBLEM * 
X ) 

RE TURN 1 
END 
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